• Uma série temporal (ST) é uma coleção de observações feitas sequencialmente no tempo.

• Muitas STs que conhecemos estão no campo da economia, as cotações diárias do mercado de ações, índice de inflação ou taxa de desemprego mensal.

• Em ciências sociais seguem séries populacionais, como taxas de natalidade.

• Na computação, a evolução de um software pode ser estudado através de suas releases.

• Na epidemiologia surge o número de casos de dengue ao longo de um período de tempo.

• Na medicina, a ressonância magnética de ondas cerebrais pode ser usada para estudar como o cérebro reage a certos estímulos experimentais.

• Na área de saúde, o comportamento de uma série temporal ficou mais conhecido com os casos de COVID-19 apresentados diariamente.

• Um dos primeiros registros de séries temporais é o número mensal de manchas solares, de 1749 a 1983, estudado por Schuster (1906).

Figura 1: Número mensal de manchas solares.

• Para reproduzir a Figura 1 no R, os comandos são:

sunspots.m <- data.frame(year = index(sunspot.month),
 value = melt(sunspot.month)$value)
ggplot(sunspots.m, aes(x = year, y = value)) +
  geom_line(colour = "red")+
  ylab("Número de manchas solares") + 
  xlab("Ano") 

• A principal diferença entre séries temporais e dados transversais está na estrutura de dependência.

• A análise de dados transversais trata, principalmente, observações i.i.d – independentes e identicamente distribuídas.

• Enquanto que em séries temporais, cada nova observação é estocasticamente dependente de observações anteriores.

• Ao mesmo tempo que a dependência é nossa maior inimiga, ela será nossa grande aliada.

• Em estudos de séries temporais é possível avaliar, através de medições de temperaturas, a hipótese de aquecimento global.

• É possível avaliar também se os níveis de poluição podem influenciar a mortalidade diária em grandes metrópolis.

• Os dois últimos são exemplos de dados experimentais que podem ser usados para ilustrar o processo pelo qual a metodologia estatística clássica é aplicada em séries temporais correlacionadas.

• De modo geral, o primeiro passo em qualquer investigação de séries temporais é a observação dos dados traçados ao longo do tempo.

• É importante mencionar que existem duas abordagens distintas, mas não mutuamente excludentes, na análise de séries temporais:

  • a abordagem no domínio do tempo;
  • a abordagem no domínio da frequência.

• A abordagem no domínio do tempo é motivada pela presunção de que a correlação entre pontos adjacentes no tempo é melhor explicada em termos de uma dependência entre o valor atual e valores passados.

• Assim, busca-se modelar valores futuros de uma série temporal como uma função paramétrica dos seus valores atuais e passados.

• Neste cenário, podemos usar regressões lineares do valor presente de uma série temporal sobre seus próprios valores passados e sobre valores passados de outras séries.

• Por outro lado, a abordagem do domínio da frequência pressupõe que as características de interesse em análises de séries temporais se relacionam com variações sistemáticas encontradas naturalmente na maioria dos dados.

• Essas variações periódicas são frequentemente causadas por fenômenos biológicos, físicos ou ambientais:

  • vento pode afetar vibrações em uma asa de avião;
  • temperaturas da superfície do mar causadas pelo El Niño podem afetar o número de peixes no oceano;

• Os estudos da periodicidade são comuns na economia e nas ciências sociais, onde há interesse na periodicidade anuais de séries como o desemprego mensal ou na taxa de natalidade mensal.

• Alguns dos problemas e questões de interesse do analista de séries temporais podem ser apresentados através de dados reais extraídos de diferentes áreas de concentração.

• Os exemplos a seguir ilustram alguns tipos de dados de séries temporais, bem como algumas questões estatísticas que podem surgir sobre esses dados.

• Esta série temporal refere-se aos dados de vendas mensais de medicamentos antidiabéticos na Austrália de 1992 a 2008.

• São 204 meses. Total de scripts mensais para produtos farmacêuticos abrangidos pelo código A10 do ATC, tal como registado pela Australian Health Insurance Commission.

• A série a10 pertence ao pacote fpp.

Figura 2: Vendas mensais de medicamentos antidiabéticos.

• Observem a tendência de crescimento e a variação regular que parece se repetir em relação aos anos.

• Para reproduzir a Figura 2 no R, os comandos são:

ggplot(a10, aes(as.Date(a10), as.matrix(a10))) + 
  geom_line(colour = "red")+
  ylab("Vendas") + 
  xlab("Tempo") 

• Esta série temporal refere-se aos dados mensais de mortes por pneumonia e gripe a cada 10.000 pessoas nos Estados Unidos por 11 anos, de 1968 a 1978.

• São 132 meses.

• A série flu pertence ao pacote astsa.

Figura 3: Mortes por pneumonia e gripe a cada 10.000 pessoas.

• Notem uma tendência leve de decrescimento e uma variação que parece se repetir em relação aos anos.

• Para reproduzir a Figura 3 no R, os comandos são:

ggplot(flu, aes(as.Date(flu), as.matrix(flu))) + 
  geom_line(colour = "red")+
  ylab("Taxa") + 
  xlab("Tempo") 

• Esta série temporal refere-se aos dados casos de AIDS notificados no SINAN, declarados no SIM e registrados no SISCEL/SICLOM(1), segundo capital de residência por ano de diagnóstico. Brasil, 1980-2021.

• Fonte: MS/SVS/Departamento de Doenças de Condições Crônicas e Infecções Sexualmente Transmissíveis (DCCI).

> AIDS <- read.csv("Dados SUS/AIDS.csv", sep = ";")
> AIDS <- ts(AIDS$Frequência, start=1980, frequency=1) 

• São 41 anos.

• A série foi coletada no site do DATASUS.

• Notem uma tendência forte de crescimento até 2013 e não apresenta variação em relação aos anos.

Figura 4: Casos de AIDS notificados no SINAN.

• Para reproduzir a Figura 4 no R, os comandos são:

ggplot(AIDS, aes(as.Date(AIDS), as.matrix(AIDS))) + 
  geom_line(colour = "red")+
  ylab("Casos") + 
  xlab("Tempo") + 
  scale_y_continuous(labels = addUnits)

• A camada scale_y_continuous(labels = addUnits) utilizada no gráfico anterior roda com base na função abaixo

addUnits <- function(n) {
  labels <- ifelse(n<1000,n,#less than thousands
             ifelse(n<1e6,paste0(round(n/1e3),'k'),#in thousands
              ifelse(n<1e9,paste0(round(n/1e6),'M'),#in millions
               ifelse(n<1e12,paste0(round(n/1e9),'B'),#in billions
                ifelse(n<1e15,paste0(round(n/1e12),'T'),#in trillions
                'too big!'
                 )))))
  return(labels)
}

• Podemos ter interesse em analisar várias séries temporais ao mesmo tempo.

• Considere as séries de um estudo para investigar possíveis associações entre mortalidade e poluição no condado de Los Angeles.

• Neste exemplo, vamos considerar as séries temporais no período de 1970 a 1979:

  • cmort Série semanal de mortalidade cardiovascular;
  • part Série semanal da quantidade de partículas de poluição.

• Ambas as séries com 508 semanas. Para maiores detalhes http://www.sungpark.net/ShumwayAzariPawitan88.pdf

• As séries cmort e part pertencem ao pacote astsa.

Figura 5: Mortalidade cardiovascular.

• Para construir as figura anteriores foram utlizados os comandos abaixo:

autoplot(cbind(Mortality=cmort), xlab='Tempo', 
         ylab='Mortalidade', col=2) 

autoplot(cbind(Mortalidade=cmort, Partículas=part), 
         xlab='Tempo', ylab='')

• Ambas as séries tendem a apresentar comportamento repetitivo, com ciclos regulares.

• Parece haver uma relação entre as duas séries.

• Esta possibilidade sugere o uso de técnicas da análise de regressão como um procedimento para relacionar as duas séries.

• Dado que temos uma ou mais séries temporais para realizar a análise, como podemos fazê-la?

• A característica peculiar de dados de séries temporais é que as observações não são independentes.

• Logo, qualquer análise deve considerar a ordem em que as observações são coletadas.

• Cada observação numa série é uma observação bivariada com o tempo sendo a segunda variável.

• Os principais objetivos da análise de séries temporais são:

1. Descrição. Descrever os dados usando estatísticas descritivas e métodos gráficos;

2. Modelagem. Encontrar um modelo estatístico adequado para descrever o processo gerador de dados. Obviamente, todos os modelos são aproximações e a construção de modelos é uma arte;

3. Previsão. Estimar (prever) valores futuros da série.

• Em análise de séries temporais é crucial desenvolver modelos estatísticos que forneçam descrições plausíveis para os dados.

• Encontrar um modelo estatístico capaz de descrever características dos dados que aparentemente flutuam de forma aleatória ao longo do tempo.

• Uma série temporal pode ser definida como uma coleção de variáveis de acordo com a ordem em que são obtidas no tempo.

• Por exemplo, uma série temporal é uma sequência de variáveis, \(x_1, x_2,\dots\), em que \(x_1\) denota o valor observado no primeiro ponto de tempo, assim por diante.

• Consideremos aqui a série temporal \(\{x_t\}\), indexada por \(t = 1, 2,\dots\), tipicamente discreto.

• Podemos imaginar que, em cada uma das séries apresentadas antes, há a suposição de que pontos adjacentes no tempo estão correlacionados.

• Ou seja, o valor no tempo \(t\), digamos \(x_t\), depende dos valores passados \(x_{t-1}, x_{t-2},\dots\)

• Para ilustrar, vamos considerar o exemplo a seguir.

• Este é um exemplo típico de série com valores não-correlacionados, \(w_t\), com média 0 e variância \(\sigma_w^2\):

\[w_t \sim wn(0, \sigma_w^2).\]

• Um caso particular é a série ruído branco gaussiano, em que \(w_t\) são independentes e normalmente distribuídas, com média 0 e variância \(\sigma_w^2\):

\[w_t \sim {\rm N}(0, \sigma_w^2).\]

• A figura a seguir mostra uma série \(wn\) gaussiano de tamanho 500, com \(\sigma_w^2 = 1\).

set.seed(1)
wn <- rnorm(500,0,1)

head(wn)

## [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684

ts_wn <- data.frame(Time = seq(0, length(wn)-1, 1), WN = wn)
head(ts_wn) #mostra as 6 primeiras observações

##   Time         WN
## 1    0 -0.6264538
## 2    1  0.1836433
## 3    2 -0.8356286
## 4    3  1.5952808
## 5    4  0.3295078
## 6    5 -0.8204684

tail(ts_wn)#mostra as 6 últimas observações

##     Time         WN
## 495  494  3.8102767
## 496  495 -1.1089100
## 497  496  0.3075666
## 498  497 -1.1068945
## 499  498  0.3476536
## 500  499 -0.8732645

• Se o comportamento estocástico de todas as séries temporais pudessem ser explicado em termos do modelo de ruído branco, os métodos estatísticos clássicos seriam suficientes.

Lembram qual o pressuposto sobre os resíduos nos modelos clássicos de regressão?

• Portanto, devemos considerar a correlação serial.

• Duas formas de introduzir a correlação serial em modelos de séries temporais são exemplificadas a seguir.

• Podemos suavisar a série de ruído branco \(w_t\) através do modelo média móvel.

• Por exemplo, considere substituir \(w_t\) no exemplo anterior por uma média de seu valor atual e seus vizinhos imediatos no passado e no futuro.

\[v_t = \frac{w_{t-1} + w_t + w_{t+1}}{3}\]

• A série \(v_t\) gerada é a da figura a seguir.

• Para reproduzir a Figura 6 e a Figura 7 no R, os comandos são:

gg_wn <- ggplot(data = ts_wn, aes(x = Time, y = WN)) + 
  geom_line() + labs(x = "\n Tempo ", y = "RB \n") +
  theme(plot.title = element_text(face="bold", 
                                  colour="darkblue", 
                                  size = 12,hjust = 0.5))
gg_wn

ts_mm <- data.frame(Time = seq(0, length(wn)-1, 1), 
                    MM = rollapply(wn,3,mean,align='right',
                                   fill=NA))
gg_mm <- ggplot(data = ts_mm, aes(x = Time, y = MM)) + 
  geom_line() + labs(x = "\n Tempo ", y = "MM \n") +
  theme(plot.title = element_text(face="bold", 
                                  colour="darkblue", 
                                  size = 12,hjust = 0.5))
gg_mm

grid.arrange(gg_wn,gg_mm, ncol = 1)

• Examinando o gráfico, vemos que \(v_t\) é uma versão mais suavizada com relação à primeira série \(w_t\).

• Ou seja, na média móvel as oscilações são mais lentas e algumas das oscilações mais rápidas são extraídas.

• Consideremos a série ruído branco \(w_t\) como entrada e vamos calcular a saída usando a equação de segunda ordem:

\[x_t = x_{t-1} - 0.9 x_{t-2} + w_t,\] para \(t=1, 2, \dots, 500\).

• A equação acima representa a regressão do valor atual \(x_t\) de uma série temporal em função dos dois últimos valores da própria série; por isso, o termo autorregressão.

• A série resultante desta autorregressão é mostrada na figura a seguir.

• Nota-se alguma periodicidade na série.

• O modelo autorregressivo é bastante utilizado no ajuste de muitas séries observadas e serão estudados em detalhes mais adiante.

• Para gerar (simular) e plotar estes dados autorregressivos no R:

set.seed(2012)
yt <- arima.sim(list(order=c(2,0,0), ar=c(1,-.9)), n=500)
ts_ar <- data.frame(Time = seq(0, length(yt)-1, 1), yt)

gg_ar <- ggplot(data = ts_ar, aes(x = Time, y = yt)) + 
  geom_line() +
  labs(x = "\n Tempo ", y = "AR \n", title = "Autoregressão \n") +
  theme(plot.title = element_text(face="bold", colour="darkblue", 
                                  size = 12,hjust = 0.5))
gg_ar

• A função de autocorrelação (ACF) é definida na forma

\[ \rho(s,t) = \frac{\gamma(s,t)}{\sqrt{\gamma(s,s) \gamma(t,t)}}. \]

• O ACF mede a previsibilidade linear da série no tempo \(t\), digamos \(x_t\), usando apenas o valor \(x_s\).

• Tem-se que \(-1 \leqslant \rho(s,t) \leqslant 1\).

• Se é possível prever perfeitamente \(x_t\) a partir de \(x_s\) através da relação linear \(x_t = \beta_0 + \beta_1 x_s\), então:

  • \(\rho(s,t) = +1\), quanto \(\beta_1 > 0\);

  • \(\rho(s,t) = -1\), quanto \(\beta_1 < 0\);

• Então, \(\rho(s,t)\) é uma medida aproximada da capacidade de prever a série no instante \(t\) com o valor no instante \(s\).

• A função de autocorrelação amostral é definida como

\[ \widehat\rho(h) = \frac{\sum\limits_{t=1}^{n-h} (x_{t+h} - \bar{x})(x_t - \bar{x})}{\sum\limits_{t=1}^{n} (x_t - \bar{x})^2}, \] para \(h = 0, 1,\dots, n-1\).

• Estimar a autocorrelação em séries temporais é semelhante a estimar a correlação no caso clássico.

• A Figura 8 mostra um exemplo usando a série de Mortalidade Cardiovascular (cmort) onde \(\widehat\rho(1) = 0.772\) e \(\widehat\rho(6) = 0.528\).

• O código a seguir foi usado para construir a Figura 8.

rho_ = round(acf(cmort, 6, plot=FALSE)$acf[-1], 3) # valores da ACF
cmort.lag1 <- dplyr::lag(as.matrix(cmort),1)
cmort.lag6 <- dplyr::lag(as.matrix(cmort),6)
df_cor <- data.frame(cmort.lag1, cmort.lag6, as.matrix(cmort))
g1<-ggplot(data = df_cor, aes(cmort.lag1, cmort)) + 
  geom_point(color='blue') + 
  geom_label(data = NULL, x = 80, y = 125, label = rho_[1], size=5, col="red") 

g2<-ggplot(data = df_cor, aes(cmort.lag6, cmort)) + 
  geom_point(color='blue') + 
  geom_label(data = NULL, x = 80, y = 125, label = rho_[6], size=5, col="red") 

grid.arrange(g1,g2, ncol=2, nrow=1)

• A função de autocorrelação amostral nos permite avaliar se os dados provêm de uma série completamente aleatória ou se as correlações são estatisticamente significantes em algumas defasagens (lags).

Propriedade: Se \(x_t\) é ruído branco, para \(n\) grande, a ACF amostral, \(\widehat\rho_x(h)\), para \(h = 1, 2,\dots, H\), é aproximadamente normal com média zero e desvio-padrão dado por \(\frac{1}{\sqrt{n}}\).

• Baseado nesta propriedade, podemos avaliar se uma série é ruído branco determinando quantos valores de \(\widehat\rho_x(h)\) estão fora do intervalo \(±1,96/\sqrt{n}\), geralmente arredondado para \(±2/\sqrt{n}\) (dois erros-padrão).

• Para o ruído branco, aproximadamente 95% dos ACFs amostrais devem estar dentro desses limites.

Importante: É devido a esta propriedade que se aplicam técnicas para transfomar uma série temporal em uma série de ruído branco, de modo a não violar a suposição de alguns modelos estatísticos.

• A ACF de cmort exibe periodicidades correspondentes a valores separados por 17 unidades de tempo, que correspondem a 17 semanas ou 1 quadrimestre.

  • Observações das semanas referentes aos \(1^{o}\) e \(3^o\)s quadrimestres de cada ano estão correlacionadas positivamente;

  • Observações das semanas referentes aos \(2^o\)s de cada ano estão negativamente correlacionadas.

• A ACF de part exibe periodicidades correspondentes a valores separados por 12 unidades de tempo, que correspondem a 12 semanas ou 1 trimestre.

  • Observações das semanas referentes aos \(1^o\)s e \(4^o\)s trimestres de cada ano estão correlacionadas positivamente;
  • Observações das semanas referentes aos \(2^o\)s e \(3^o\)s trimestres de cada ano estão negativamente correlacionadas.

• Para obter estas ACFs das Figuras Figura 9 e Figura 10 no R, usem os comandos:

ggAcf(cmort, 104, main="Mortalidade Cardiovascular")
ggAcf(part, 104, main="Partículas de Poluição")

• Até aqui, não fizemos qualquer suposição especial sobre o comportamento da série temporal.

• Os exemplos anteriores sugerem alguma regularidade ao longo do tempo no comportamento da série.

• Assim como nos modelos clássicos de regressão, quando pretendemos utilizar modelos para descrever séries temporais, é necessário introduzir suposições simplificadoras.

• A principal suposição é a de estacionariedade.

Intuitivamente, uma série \(x_t\) é estacionária se ela se desenvolve no tempo de modo que a escolha de um instante no tempo não é importante (Morettin e Toloi, 2006).

• Em outras palavras, as características de \(x_{t+s}\), para todo \(s\), são as mesmas de \(x_t\).

• Tecnicamente, há duas formas de estacionariedade:

  • fraca (ou de segunda ordem); e
  • estrita (ou forte).

• Para mais detalhes, sugiro ver Morettin e Toloi (2006).

• O modelo de regressão linear é mais utilizado no contexto de dados clássicos do que em séries temporais.

• Estes modelos de regressão são importantes para compreender os modelos de séries temporais que discutiremos mais adiante.

• Já vimos que os modelos autorregressivos expressam \(x_t\) como uma combinação linear de valores passados \(x_{t-1}, x_{t-2},\dots, x_{t-p}\) da própria série.

• Além disso, \(x_t\) pode depender de valores defasados de outra série, digamos \(y_{t-1}, y_{t-2},\dots, y_{t-q}\), que a influencia.

• Assumimos que a série temporal dependente \(x_t\), para \(t = 1,\dots, n\), é influenciada por um conjunto de séries independentes, \(z_{t1}, z_{t2},\dots, z_{tq}\), conhecidas.

• Esta é uma condição necessária para aplicar a regressão linear convencional.

• Então, expressamos essa relação através do modelo de regressão linear

\[ x_t = \beta_1 z_{t1} + \beta_2 z_{t2} + \cdots + \beta_q z_{tq} + w_t, \] em que \(\beta_1, \beta_2,\dots, \beta_q\) são os coeficientes desconhecidos a serem estimados, \(w_t\) é o erro aleatório (ou ruído) independente e normalmente distribuído com média 0 e variância \(\sigma_w^2\).

• Voltemos ao exemplo dos Casos de AIDS na figura a seguir.

• Lembrem que nesta série há uma aparente tendência ascendente que é usada para argumentar a hipótese do aumento dos casos de AIDS.

• Podemos usar regressão linear para estimar essa tendência ajustando o modelo simples

\[ x_t = \beta_1 + \beta_2 t + w_t,\qquad t = 1980, \dots, 2020. \]

• Este é um modelo de regressão, onde \(t\) é a variável independente.

• Observe que estamos fazendo a suposição de que os erros, \(w_t\), são uma sequência iid normal.

# regressão dos casos de AIDS no tempo
modelo = lm(AIDS~time(AIDS))
summary(modelo)$coefficients

##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) -2205365.209 233235.3452 -9.455536 1.218968e-11
## time(AIDS)      1115.431    116.6156  9.565020 8.878637e-12

• Pela regressão linear simples, o modelo estimado é

\[\widehat{x}_t = -2205365.21 + 1115.431\ t.\]

• O modelo estima um aumento significativo de 1115.4 casos de AIDS a cada ano.

• A figura abaixo mostra a série e a tendência estimada.

Figura 11: Casos de AIDS com tendência linear estimada.

• Observem que a linha de tendência estimada via regressão linear não capta bem a tendência dos dados.

• Abaixo temos os comandos usados para contruir a Figura 11.

ggplot(data = AIDS, aes(as.Date(AIDS), as.matrix(AIDS))) + 
  geom_point(color='blue') +
  labs(x = "\n Tempo ", y = "Casos \n") +
  scale_y_continuous(labels = addUnits)+
  geom_smooth(method = "lm", se = FALSE, color='red')

• Em geral, é desejável que a série temporal seja estacionária, de modo que a sua média não mude com o tempo.

• Além disso, a dependência entre os valores da série é algo que devemos medir.

• No entanto, é difícil medir esta dependência se a estrutura dela estiver mudando com o tempo.

• Então, para obter qualquer análise estatística significativa da série, é crucial que a média e a função de autocovariância satisfaçam as condições de estacionaridade:

  1. a média, \(\mu_t\), é constante e não depende do tempo \(t\), e

  2. a função de autocovariância, \(\gamma(s,t)\), depende de \(s\) e \(t\) somente através de sua diferença \(\vert s-t \vert\).

• Muitas vezes, este não é o caso, e vamos mencionar alguns métodos nesta seção para minimizar os efeitos da não-estacionaridade.

• A série de casos de AIDS contém evidência de uma tendência ao longo do tempo:

• Uma forma comum de lidar com a não-estacionaridade é quando o processo tem comportamento estacionário em torno de uma tendência.

• Podemos escrever este tipo de modelo como

\[x_t = \mu_t + y_t,\] em que \(x_t\) é a série, \(\mu_t\) denota a tendência e \(y_t\) é um processo estacionário.

• Logo, pode ser vantajoso remover a tendência como um primeiro passo em uma análise exploratória de séries temporais.

• Para isso, basta obter uma estimativa razoável do componente de tendência, digamos \(\widehat{\mu}_t\), e então trabalhar com os resíduos

\[ \widehat{y}_t = x_t -\widehat{\mu}_t. \]

• A tendência estimada via MQO para esta a série das vendas de medicamentos antidiabéticos foi

\[ \widehat{\mu}_t = -2205365.21 + 1115.431\ t. \]

• Para obter a série sem tendência, nós simplesmente subtraímos \(\widehat{\mu}_t\) das observações \(x_t\)

\[ \begin{aligned} \widehat{y}_t &= x_t -\widehat{\mu}_t\\[1mm] \widehat{y}_t &= x_t + 2205365.21 -1115.431\ t. \end{aligned} \]

• A nova série \(\widehat{y}_t\) com tendência linear removida é apresentada na figura seguinte.

• Abaixo temos os comandos usados para contruir a Figura 12.

modelo = lm(AIDS~time(AIDS), na.action=NULL)
residuo = resid(modelo)
ggplot(data = residuo, aes(as.Date(residuo), as.matrix(residuo))) + 
  geom_line(color='blue') +
  labs(x = "\n Data ", y = "Resíduos \n") 

• A figura abaixo mostra a ACF amostral para a série original dos casos de AIDS.
  

• E a ACF amostral para a série dos casos de AIDS com tendência removida.

• Para eliminar a tendência no R, plotar a série sem tendência e as ACF amostrais anteriores, usamos os seguintes comandos:

modelo = lm(AIDS ~ time(AIDS), na.action=NULL)
ggplot(data = AIDS, aes(as.Date(AIDS), as.matrix(AIDS))) + 
  geom_point(color='blue') +
  labs(x = "\n Data ", y = "Casos \n") 
plot(resid(modelo), type="o", main="sem tendência")
par(mfrow=c(3,1)) # plot ACFs
ggAcf(AIDS, 48, main="Casos de AIDS")
ggAcf(resid(modelo), 48, main="sem tendência")

• Outra prática comum para tornar uma série estacionária é tomar as diferenças sucessivas.

• Uma vantagem de diferenciar a série sobre o processo de remover tendência anterior é que nenhum parâmetro é estimado na operação de diferenciação.

• A primeira diferença é denotada por

\[ \nabla x_t = x_t - x_{t-1}. \]

• Em geral, a primeira diferença elimina uma tendência linear.

• Uma segunda diferença, dada por

\[ \begin{aligned} \nabla^2 x_t &= \nabla(\nabla x_t)\\[1mm] &= \nabla(x_t - x_{t-1})\\[1mm] &= (x_t - x_{t-1}) - (x_{t-1} - x_{t-2})\\[1mm] &= x_t - 2x_{t-1} + x_{t-2}, \end{aligned} \] pode eliminar uma tendência quadrática, e assim por diante.

• De modo geral, a \(d\)-ésima diferença de \(x_t\) é

\[ \nabla^d x_t = \nabla(\nabla^{d-1} x_t). \]

• Em situações normais, será suficiente tomar uma ou duas diferenças para que a série se torne estacionária.

• A técnica de diferenciação é um componente importante do modelo ARIMA discutido mais adiante.

• A primeira diferença da série dos casos de AIDS produz resultados diferentes da remoção de tendência via regressão.

autoplot(diff(AIDS), type="o", main="primeira diferença")

• A ACF amostral da série dos casos de AIDS diferenciada é

• Parece que a série da primeira diferença já mostra autocorrelação mínima.

• Para obter e plotar os dados da primeira diferença, e gerar a ACF amostral, usamos os comandos:

plot(diff(AIDS), type="o", main="primeira diferença")
ggAcf(diff(AIDS), 48, main="primeira diferença")

• Esta é outra técnica exploratória de dados usada para visualizar relações entre diferentes defasagens de séries temporais.

• Anteriormente, vimos que a ACF nos indica existência de relação linear entre a série e seus próprios valores defasados.

• No entanto, a ACF pode mascarar uma possível relação não-linear entre \(x_t\) e seu valores passados, \(x_{t-h}\).

• Para avaliar possíveis relações não-lineares via matriz de dispersão, podemos plotar a série original, \(x_t\), contra suas defasagens \(x_{t-h}\).

# install.packages("astsa")
library(astsa)
lag1.plot(cmort, 12)

• As autocorrelações amostrais são exibidas no canto superior direito nos gráficos de dispersão.

• Os gráficos também exibem as linhas de alisamento de dispersão (lowess) usadas para indicar a relação não-linear.

Observação: Discutiremos alisamento a seguir, mas entenda lowess como um método para ajustar a regressão local.

• Na figura anterior, observamos que os ajustes lowess são aproximadamente lineares, com autocorrelações significativas.

• Além disso, vemos fortes relações lineares positivas entre \(M_t\) e as defasagens \(M_{t-1}, M_{t-2}, M_{t-3}, M_{t-4}\) e depois começa a reduzir o valor da correlação à medida que aumenta a quantidade de defasagens.

• Da mesma forma, podemos avaliar possíveis relações não-lineares entre a série de mortalidade cardiovascular, \(M_t\), e desfasagens do nº de partículas, \(P_{t-h}\).

Lembrem do exemplo onde queríamos prever a média de mortalidade semanal a partir da quantidade de partículas de poluição.

• A figura anterior mostra a matrix de dispersão da série de mortalidade cardiovascular contra a série do nº de partículas e suas defasagens, \(P_{t-h}\).

• Além disso, a figura exibe as correlações cruzadas, bem como os ajustes lowess.

• Nota-se relação não-linear de forma moderada entre mortalidade cardiovascular \(M_t\) e partículas de poluição nas defasasgens \(P_{t-4}, P_{t-6}, P_{t-7}, P_{t-8}\), indicando que a série de partículas pode influenciar antecipadamente a série de mortalidade cardiovascular.

• Para construir as figuras das matrizes de dispersão, usamos o pacote astsa e o código no R é

install.packages("astsa")
library(astsa)
astsa::lag1.plot(part, 12)        # dispersão com uma série
astsa::lag2.plot(part, cmort, 8) # dispersão com duas séries

• No Exemplo: Média móvel discutimos a utilização de média móvel para suavizar o ruído branco.

• Este método é útil para descobrir algumas características em séries temporais, como tendências de longo prazo e componentes sazonais.

• Em particular, se \(x_t\) representa as observações, então \[ m_t = \sum\limits_{j=-k}^{k} a_j x_{t-j}, \] em que \(a_j = a_{-j} \geqslant 0\) e \(\sum_{j=-k}^{k} a_j = 1\) é a média móvel dos dados.

• Para ilustrar, a figura a seguir mostra a série de Óbitos por COVID-19 na Bahia, suavizada via média móvel

Figura 15: Óbitos diários por COVID-19 na Bahia

• Comandos usados:

#library(readxl)
covid_bahia <- 
  read_excel("Dados/Banco Estadual COVID-19 ÓBITOS_17-10-2022.xlsx")

covid_bahia$Mortes <- seq(1,1,length(covid_bahia$`DATA OBITO`))
#View(covid_bahia)

class(covid_bahia$`DATA OBITO`)

• Comandos usados:

### Ajustando a variável data usando 
##  o pacote lubridate
## Indicar a ordem da data em inglês 
## Y-year , m- month d-day
#library(lubridate)
covid_bahia$`DATA OBITO` <- 
  as.Date(parse_date_time(covid_bahia$`DATA OBITO`, 
                          "%y/%m/%d"))
class(covid_bahia$`DATA OBITO`)
df_covid = tibble(covid_bahia)

dplyr::glimpse(df_covid)

• Comandos usados:

df_covid <- df_covid%>%
  dplyr::group_by(`DATA OBITO`)%>%
  dplyr::summarise_if(is.numeric,sum)

#library(zoo)
df_covid=df_covid %>%
  mutate('Média Móvel'=rollapply(Mortes,14,
                                 mean,align='right',
                                 fill=NA))

• Comandos usados:

## Agrupando para termos categorias de 
## Médias Móvel e Óbitos no mesmo gráfico
### Deixando em outro conjunto Bahia

df_covid2=df_covid %>%
  gather(c('Média Móvel','Mortes'),key="Séries",
         value="Valor")

ggplot(df_covid2,aes(x=`DATA OBITO`,y=Valor,fill=Séries, 
                     colour=Séries))+
  geom_line()+
  labs(x="",y="Número de óbitos")

• Outra abordagem para suavizar uma série é a regressão do vizinho mais próximo.

• A técnica baseia-se na regressão de \(k\)-vizinhos mais próximos, onde se usa apenas os dados \(\{x_{t-k/2},\dots, x_t,\dots, x_{t+k/2}\}\) para prever \(x_t\) via regressão, e então define \(m_t = \widehat{x}_t\).

Lowess é um método de alisamento complexo, mas a ideia básica é próxima à regressão do vizinho mais próximo.

• A figura a seguir mostra a suavização da série de óbitos por COVID usando a função lowess no R (ver Cleveland, 1979).

• Código da figura anterior.

• A camada geom_smooth() tem como default o método lowess

ggplot(df_covid,aes(x=`DATA OBITO`,y=Mortes))+
  geom_point(col="darkblue") +
  geom_smooth(col="red")

• Podemos usar o argumento span para controlar a suavização.

• Números menores produzem linhas mais onduladas, números maiores produzem linhas mais suaves.

ggplot(df_covid,aes(x=`DATA OBITO`,y=Mortes))+
  geom_point( col="darkblue") +
  geom_smooth(col="red")+
  geom_smooth(span=0.08, col="yellow", size=1.1)

• Comandos usados para construir a figura anterior.

• Usando o exemplo da série com os casos de AIDS.

• Comandos usados para construir a figura anterior.

ggplot(AIDS,aes(x=time(AIDS),y=as.matrix(AIDS)))+
  geom_point(col="darkblue") +
  geom_smooth(col="red")

• Técnicas de suavização podem ser aplicadas também no alisamento de uma série temporal em função de outra.

• Aqui, vamos suavizar o diagrama de dispersão de duas séries temporais medidas simultaneamente.

• Vamos considerar os dados de possíveis efeitos da temperatura (temp) sobre a mortalidade cardiovascular (cmort) em Los Angeles.

• A figura seguinte mostra o diagrama de dispersão de mortalidade, \(M_t\), vs. temperatura, \(T_t\), e a linha suavizada via lowess.

• Notem uma relação não-linear entre mortalidade e temperatura.

• A mortalidade aumenta a temperaturas extremas, mas de uma forma assimétrica; é mais alta em temperaturas mais frias do que em temperaturas mais quentes.

• A taxa de mortalidade mínima ocorre a aproximadamente 83\(^{\circ}\) F.

• O código usado para reproduzir esta figura no R foi:

df_tm <- data.frame(tempr,cmort)

ggplot(data = df_tm, aes(tempr, cmort)) + 
  geom_point(color='blue') + 
  labs(x = "\n Temperatura (em Fº) ", y = "Mortalidade \n") +
  geom_smooth(col="red")

• Decomposição de séries é um procedimento matemático útil para dividir (decompor) a série em três componentes:

  • Tendência;

  • Sazonalidade;

  • Ruído (resíduo) aleatório.

• Há dois tipos clássicos de decomposição; a aditiva,

série temporal = tendência + sazonalidade + ruído aleatório

e a multiplicativa,

série temporal = tendência * sazonalidade * ruído aleatório.

• Para obter uma boa decomposição da série, o primeiro passo é escolher entre o modelo aditivo ou multiplicativo.

• Duas regras básicas:

  • Se a variação sazonal parece constante ou não varia muito, aplicamos o modelo aditivo;

  • Se a magnitude da sazonalidade aumenta/diminui com o tempo, devemos aplicar o modelo multiplicativo.

• Para ilustrar a técnica de decomposição, consideremos as duas séries temporais a seguir.

• A primeira é a série da produção trimestral de cerveja na Austrália (em megalitros); 70 observações a partir do 1º trimestre de 1956.

Figura 17: Produção de cerveja na Austrália.

Série da produção de cerveja na Austrália

• A variação sazonal não varia com o tempo, i.e., é constante.

• Então, vamos usar o modelo aditivo para esta série.

• A segunda é a série de vendas mensais de medicamentos antidiabéticos na Austrália de 1992 a 2008.

Série de vendas mensais de medicamentos antidiabéticos

• A variação sazonal aumenta ao longo do tempo.

• Logo, devemos usar o modelo multiplicativo para esta série.

• Uma função do R útil para visualizar a variação sazonal de uma série é seasonplot do pacote forecast.

seasonplot(st.ausbeer, s=4, col=rainbow(4))
seasonplot(a10, s=12, col=rainbow(12))

• Ela produz o gráfico sazonal, útil na escolha do modelo a ser usado na decomposição.

• Este gráfico sazonal permite visualizar o padrão sazonal da série a cada ciclo, e.g., a cada ano, e é útil na identificação de anos em que o padrão muda.

• Observem, por exemplo, o aumento na produção de cervejas na Austrália no 2º trimestre de cada ano.

• Para as vendas de medicamentos, notem o aumento, principalmente no mês de janeiro.

• Para ilustrar a técnica de decomposição, vamos usar algumas funções do pacote TSstudio e reproduzí-la passo-a-passo para entender como ela funciona.

• Faremos o exemplo para demonstrar o passo-a-passo com a série que apresenta sazonalidade aditiva, st.ausbeer.

[Passo 1]

• Inicialmente, vamos detectar a tendência das séries exemplificadas, via suavização média móvel.

• Para decompor uma série, é crucial conhecer o período de sazonalidade: semanal, mensal, anual, etc.

• A produção de cerveja australiana tem sazonalidade anual, ou seja, os ciclos parecem se repetir a cada ano.

• Como é registrada trimestralmente, há 4 observações por ano, então, usamos janelas (vizinhos imediatos) de média móvel de 4.

[Passo 1]

• A regra geral é usar média móvel de dois lados com ordem de frequência/2.

• Por exemplo, no caso do dataset st.ausbeer, como a frequência é trimestral, vamos usar média móvel de dois lados, onde cada lado terá 2 observações consecutivas.

• Isso significa que, para calcular o valor suavizado da observação \(t\), calcularemos a média da observação \(t\) junto com as 2 observações anteriores e seguintes (ou seja, com média de 5 pontos dos dados).

• Estimamos a tendência com a função ts_ma {TSstudio} do R:

ap_smooth <- ts_ma(st.ausbeer, n = 2, plot = FALSE, 
                   separate = FALSE)
ap_smooth$plot %>%
  layout(legend = list(x = 0.1, y = 0.9))

• Na próxima etapa, estimamos a tendência

• Converteremos a série e a tendência suavizada em um objeto data.frame com a função ts_to_prophet para mesclar as duas séries:

# modelo aditivo
ap_smooth <- ts_ma(st.ausbeer, n = 2, plot = FALSE, separate = FALSE)
df <- ts_to_prophet(st.ausbeer) %>% 
  select(date = ds, y) %>% 
  left_join(ts_to_prophet(ap_smooth$ma_2) %>%
              select(date = ds, trend = y), by = "date")

df$detrend <- df$y - df$trend
df$year <- year(df$date)
df$month <- month(df$date)

ts_plot(df[, c("date", "y", "trend")], title = "Beer production and trend component",
        type = "multiple")

[Passo 2]

• Agora, vamos remover a tendência previamente calculada resultando em uma nova série temporal que exibe claramente a sazonalidade.

# modelo aditivo

df$detrend <- df$y - df$trend

df$year <- year(df$date)
df$month <- month(df$date)
p <- plot_ly()
for(i in unique(df$year)){
  temp <- NULL
  temp <- df %>% filter(year == i) 
  p <- p %>% add_lines(x = temp$month,
                       y = temp$detrend,
                       name = i)
  
}

seasonal_comp <- df %>% 
  group_by(month) %>%
  summarise(month_avg = mean(detrend, na.rm = TRUE),
            .groups = "drop")

p %>% add_lines(x = seasonal_comp$month, 
                y = seasonal_comp$month_avg,
                line = list(color = "black", dash = "dash", width = 4),
                name = "Seasonal Component")

[Passo 3]

• Vamos, então, obter a sazonalidade média da série.

• Observem que em cada ano há um ciclo diferente na série das sazonalidades.

• Cria-se uma matriz em que as linhas representam os anos e as colunas contêm os elementos do mesmo período (mesmo mês, mesmo trimestre, etc) e calcula a média de cada coluna.

• Este passo 3 corresponde a traçar a curva sazonal média.

• Comandos

# modelo aditivo
seasonal_comp <- df %>% 
  group_by(month) %>%
  summarise(month_avg = mean(detrend, na.rm = TRUE),
            .groups = "drop")

df <- df %>% left_join(seasonal_comp, by = "month")

df$irregular <- df$y - df$trend - df$month_avg
ts_plot(df[, c("date", "month_avg")], 
        title = "Seasonal Component",
        type = "multiple")

• Graficamente

[Passo 4] Computando o ruído aleatório.

As etapas anteriores já extraíram muitas informações da série temporal original. Resta apenas o que chamamos de ruído aleatório.

Pela formulação do modelo, sabemos que na decomposição aditiva,

ruído aleatório = série temporal - tendência - sazonalidade 

e na multiplicativa,

ruído aleatório = série temporal / (tendência * sazonalidade).

Ou seja, o ruído aleatório é o que resta de informação após eliminação da sazonalidade e da tendência da série temporal.

df <- df %>% left_join(seasonal_comp, by = "month")


df$irregular <- df$y - df$trend - df$month_avg
head(df)

• Agora podemos ver a série e todas suas componentes através dos comandos:

ts_plot(df[, c("date", "y" ,"trend", "detrend", "month_avg", "irregular")], 
        title = "Beer Production and its Components",
        type = "multiple")

• Para tornar nossa vida mais fácil, alguns pacotes do R fornecem a decomposição com uma única linha de comando.

• Por exemplo decompose, stl, STL.

• As etapas apresentadas aqui reproduzem a função ts_decompose do pacote TSstudio.

• Basta indicar o período correspondente ao ciclo sazonal da série e o modelo a ser aplicado.

• Vejamos:

ts_decompose(st.ausbeer, type = "additive")
ts_decompose(a10, type = "multiplcative")

ts_covid <- stats::ts(df_covid$Mortes, start=2020-03-28, frequency=365) 
ts_decompose(ts_covid, type = "both")

• É uma medida de quão “fácil” a série é de prever.

• De acordo com Hyndman e Athanasopoulos (2021), uma série com alta sazonalidade e tendência (bem comportada) terá uma entropia perto de 0, enquanto uma série com muitos ruídos, logo difícil de prever, terá uma entropia perto de 1.

• Para obter essa medida basta usar a função feat_spectral do pacote feasts.

feat_spectral(st.ausbeer)

## spectral_entropy 
##        0.3754978

feat_spectral(a10)

## spectral_entropy 
##        0.3499013

feat_spectral(ts_covid)

## spectral_entropy 
##        0.4194675

feat_spectral(cmort)

## spectral_entropy 
##        0.7684916

• Pacotes: tidyverse, ggplot2, gridExtra, lubridate, dplyr, reshape2, zoo, readxl, astsa, fpp, forecast, magrittr, rlang, tidyverse, lubridate, feasts, zoo, fpp, TSstudio, plotly.

Acesse: https://github.com/rladiessalvador

Siga: @rladiessalvador